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Binding energy of the positronium negative ion via dimensional scaling 

Nikita Blino\Q and Andrzej Czarnccki 
Department of Physics, University of Alberta, Edmonton, AB T6G 2E1, Canada 

We determine the binding energy of the negative positronium ion in the limits of one spatial 
dimension and of infinitely many dimensions. The numerical result for the one-dimensional ground 
state energy seems to be a rational number, suggesting the existence of an analytical solution for 
the wave function. We construct a perturbation expansion around the infinitely-dimensional limit 
to compute an accurate estimate for the physical three-dimensional case. That result for the energy 
agrees to five significant figures with variational studies. 
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I. INTRODUCTION 

The negative positronium ion Ps^ is a bound state of two electrons and a positron. It is the simplest bound three 
body system from the theoretical point of view, since it does not contain a hadronic nucleus. It provides an important 
testing ground for quantum electrodynamics (QED), which should be able to describe this purely leptonic bound state 
with high precision. 

Because of the e^e~ annihilation, Ps~ is unstable, with a lifetime of about four times that of para-positronium. It 
decays predominantly into two or three photons, with the one-photon decay possible but extremely rare. It is weakly 
bound and has no excited states in the discrete spectrum 1, 2] (for a discussion of resonances, see 0, 

Following its prediction by Wheeler in 1946 |5j] and experimental observation in 1981 by Mills @, the positronium 
ion has been subject to much theoretical study. Its non-relativistic bound state energy, decay rate, branching ratios 
of various decay channels, and polarizabilities have been computed accurately using variational methods |7Hl4l|. 

Recently, intense positronium sources have become available, opening new possibilities for experimental studies of 
Ps~ The measured decay rate [IB, ES] agrees with the theoretical prediction. Improved measurements of the 

decay rate, the three-photon branching ratio, and the binding energy have been proposed fl85. 

A challenge in the theoretical study of this three-body system is that its wave function is not known analytically, 
even if only the Coulomb interaction is considered. Since all particle masses and magnitudes of their charges are 
equal, it is not possible to use the Born-Oppenheimer approximation. So far all precise theoretical predictions of Ps~ 
properties have relied on variational calculations. 

In the present paper we explore a different approach to computing the wave function and the binding energy of 
Ps~. We use dimensional scaling (DS) method, in which the dimensionality of space D is a variable. We focus on 
the limits D ^ 1 and Z? — )■ oo. A precise result for £) = 3 may be obtained by interpolating between the two limits 
using perturbation theory in l/D. The advantage of DS is that the two limits of the Schrodinger equation often 
have relatively simple solutions. Full inter-particle correlation effects are included at every order in the perturbation 
expansion in l/D. More information about dimensional scaling and further references can be found in [19l - t24| . 

It is important to note that the dimensional limits considered here are not physical in the sense that the form 
of the potential energy is taken to be l/r, regardless of the dimension. A physical limit of a system would use an 
appropriate Coulomb potential that is the solution of a D-dimensional Poisson equation. For example for D = 1 it 
is linear, logarithmic for D — 2, and depends on charge separation as r~^^~^^ for D > 2. Since we are ultimately 
interested in D = 3 physics, it is useful to fix the potential to be the D = 3 Coulomb interaction. The D — )■ 1 
limit used here offers the additional simplification that after coordinate and energy rescaling the potential takes form 
{D — l)/r , which can be formally replaced by a Dirac delta function p5| . 

We find that the DS provides a useful complement to the variational method. In the future, it can be employed to 
independently check matrix elements of operators needed in precise studies of Ps~. 

This paper is organized as follows. In Section [ll] we consider the D = 1 limit of the Ps^ system. We solve the 
Schrodinger equation numerically to find an eigenvalue that approaches a simple rational number, possibly hinting at 
the existence of an analytical solution. 

In Section Hn] we consider the _D — )■ oo limit and describe the resulting l/D expansion. We sum up the perturbation 
series for the ground state energy and evaluate it at Z? = 3. The binding energy we find agrees with variational studies 
to five significant figures. We conclude in Section ITVl 
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II. D = l LIMIT OF Ps" 



In the one dimensional limit, the Coulomb potential is represented by the Dirac delta function [25[. Delta function 
models have been used extensively also in condensed matter physics. A simple analytical wave function exists for any 
number of identical particles interacting via attractive potentials (26j . The case of all repulsive potentials with periodic 
boundary conditions has been treated by Lieb and Liniger [23|, and Yang [2^. More recent works have studied one 
dimensional systems with both attractive and repulsive delta interactions. Craig et al. considered the dependence of 
the energy on the number of particles in a system of equal numbers of positively and negatively charged bosons [2^ . 
Li and Ma studied a system of N identical particles with an impurity with periodic boundary conditions [s^l • 

The D = 1 limit of the Ps" quantum problem is a delta function model with two attractive and one repulsive delta 
functions with non-periodic boundary conditions, which, to the best of our knowledge, has not yet been solved. We 
present a derivation of a one dimensional integral equation for the solution to this problem, analogous to the helium 
case treated by Rosenthal 

The time independent Schrodinger equation for the relative motion of Ps^ takes the dimensionless form 

4 [V? + V2 + Vi.V2] ---- + -) V = £^, (1) 
2 ri r2 r J 

where ri and r2 are the electron-positron distances, r is the inter-electron distance (in units of 2{ma)^^ with h — c ~ \) 
and e determines the energy eigenvalue, E = ema^/2. This choice of units helps compare intermediate results with 
Rosenthal's delta function model of helium [sij. 

In the limit Z? — )• 1, we let ri — > x and ^2 — >■ y, where —00 < x, y < 00; the gradients become partial derivatives 
and the Coulomb potentials are replaced by Dirac delta functions (this limit is described in detail in [12]). Equation 
([T]) is replaced by 



-■1 



S{x) - S{y) + S{x - y) 



2 \ dx^ dy"^ dxdy J 

= e^. (2) 

Using Fourier transformation, we rewrite this Schrodinger equation as a one-dimensional integral equation, 

^ \{kl^kl + k^k.+p^) ' 
where the Fourier transforms of the wave function "0(3:, y) are 

G{k^M) = JJ e-''''^-'''-y^{x,y)dxdy, (4) 

F{k) = J e-'^'-'ijj{x,Q)dx, (5) 

H{k) = j e~'^'-'ip{x,x)dx, (6) 

and p2/2 = — e. We now invert the transformation and use the resulting -0 in eqs. (|5I6I) to obtain a system of two 
integral equations for F{k) and H{k). These are easily decoupled and yield 

^(fc)^ 2F(fc) , 1 [ F{k')dk' 



v/3fc2 + 4p2 n J k^ + k'^ + kk' + p"^ 
2 f ^3fc'2 + 4p2 1 / r F{k")dk 



■K^ J 2 + A/3fc'2 + 4p2 k^ + k'^ - kk' + p^ \J k'^ + fc"^ - k'k" + p^ 

and 



dk', (7) 



^ 2 v/3fc2+4p2 f F{k')dk' 

^' TT 2 + V3fc2 + 4p2 7 k^ + k'-^-kk'+p^- ^' 

Once F{k) is found one can compute H{k). The two-dimensional eigenvalue problem is thus reduced to a one- 
dimensional integral equation ([7]), which we solve numerically. The integral equation is discretized using Gauss- 
Legendre quadrature, casting it into a system of homogeneous linear equations for F{ki), where ki are the abscissas. 
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Figure 1. (Color online) The unnormalized wave function xp{x,y) satisfying eq. ((2)1. The distances are in units 2/{rna). 



The system has a non-trivial solution vifhen the determinant of the discretized integral kernel vanishes. This condition 
fixes the value of p and thus the D = 1 binding energy. 

The wave function is then determined by solving the linear system for F{ki). One finds that F{ki) spans the nuh 
space of the discretized kernel and can be computed using its singular value decomposition. We used cubic spline 
interpolation on the set {F{ki)} to interpolate between the quadrature points and generate an approximation for 

m- 

Once F{k) is known, the functions H{k) and G{ki,k2) are constructed using eq. ([5]) and Finally, the wave 
function ip{x,y) is obtained by the inverse Fourier transformation of G{ki,k2)- 

We performed this procedure for various quadrature sizes N with the results summarized in Table HI 



Quadrature size A'^ e 

10 -0.6666657902370426 

20 -0.6666666661283767 

50 -0.6666666666666257 

100 -0.6666666666666660 



Table I. Binding energy of the one-dimensional model of the positronium ion, in units mQ^/2. For N > 100 the eigenvalue 
appears to converge to —2/3. For these large quadrature sizes the uncertainty in the energy is in the last digit due to finite 
precision used in the calculation. 

We see that as N increases, e approaches —2/3. For N — 100 the 16 decimal place precision limit of the double 
data type used in the calculation is almost reached. This simple numerical result suggests that the one dimensional 
Schrodinger equation has an analytical solution. The wave function and its Fourier transform are plotted in Figures 
[T]and[2] We observe that the wavefunction has ridges at x = 0, y = and x = y a.s expected from the delta function 
potential of eq. A simple numerical comparison of H{k) with the Fourier transform of exp(— a|a;|) indicates that 
the wavefunction fall-off in the x = y direction is nearly exponential. 

The result e = —2/3 translates into the energy eigenvalue E equal —1/3 atomic unit of energy (1 a.u. = ma^) or 
—9.07 eV. This is in qualitative agreement with the actual value that is about —0.26 a.u., just below —1/4 a.u. (this 
fraction is the binding energy of a positronium atom in the non-relativistic approximation) . 

We note that the eigenvalue that can be obtained for the two-body problem (the positronium or the hydrogen 
atom) in the one dimensional delta model coincides precisely with the physical value. This is the case because the 
wave function in the delta model has the same cusp at the origin as the radial wave function in the physical space. 
Thus the delta model reproduces that radial wave function exactly. For the three-body problem the agreement is only 
rough. 

It would be interesting to determine the one dimensional wave function analytically. We remark that the Schrodinger 
equation of the three body problem ^ can be rewritten, with a simple change of variables x, y, in the form of a 
one-particle motion in the external potential consisting of two attractive and one repulsive delta function ridges. 

In the following section we focus on the opposite limit of very many dimensions. We shall find that an expansion 
around that limit can be constructed, giving a very accurate determination of the binding energy of Ps~. Interestingly, 
the D = I method will be again useful: it will provide an important subtraction term that we will use to accelerate 
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Figure 2. (Color online) The Fourier transform G{ki,k2) of the ground state position-space wave function xp{x,y) of eq. ([2]) 
computed numerically. The Fourier space coordinates have units ma/2. 



the convergence of a perturbative expansion. 



III. D ^oo AND DIMENSIONAL PERTURBATION THEORY 



The first step in taking the — >■ oo Umit is to generalize the Ps Schrodinger equation to D dimensions. We are 



interested in the ground state, which is completely described by the three inter-particle distances pij = \ri 
Schrodinger equation takes the form 

H4,= {T + U + V)(l) ^ E(j), 

where E is the energy in atomic units ma^ (note that it differs by a factor 1/2 from the e used in D 
previous section) and 



r,|. The 
(9) 

1 in the 



T 



^ I dpi 



'ipi]Pik dpijdpi 



U 
V 



P-l)(g-5) , 2 

1 



P23 + P12) 



and 



1 1 1 

P12 P13 P23 ' 

is the rescaled wave function with 
T 



s{s - pi2){s - pi3)(s - P23), 



S = 2(^^12 + Pl3 + P23) 



(10) 

(11) 
(12) 



Note that the characteristic dimensional dependence is confined to U . (The U term in the effective potential is 
the usual centrifugal contribution from the kinetic energy found by expressing the Laplacian in terms of pij.) In order 
to obtain a finite limit the coordinates and the energy must be rescaled, pij = D^Vij and E = e/D^. This introduces 
a factor of in front of the kinetic energy term, eq. pH)) . so in the limit 13 —> 00 it is suppressed. In terms of the 

rescaled quantities, the Schrodinger equation is written as 



{5^T + 5^U + V)(t> = e(j>, 



(13) 



where 5 — 1/D. 

In the limit 6^0, terms containing derivatives vanish in eq. (I13p . Since the ground state energy is the smallest 
eigenvalue of the Hamiltonian, we seek to minimize the effective potential 



VeS = S^U + V 



(14) 



5 



at (5 = 0, under the constraint that define a triangle. Unfortunately, in D — >■ oo, the Ps~ system described by 
the potential in eq. (|10p is unbound (even if the positron were very heavy, its charge would have to be larger than 
1.228 for a bound state to exists [1^ (see also [33l. Is^)). The qualitative explanation of this is that even though we 
have increased the number of spatial dimensions, we have retained the l/r behavior of the Coulomb potential. Thus 
it is relatively stronger at large distances than in three dimensions and the electron-electron repulsion plays a more 
important role even if the electrons are on the opposite sides of the positron. 

However, the strict (5 = regime is unphysical. We are interested in the 5=1/3 case, so we are free to modify 
the potential as long as it reduces to the correct form at 6 = 1/3. This can be done by reducing the strength of the 
electron-electron repulsion, as was done for H~ [23j . 

^^_ Xo + 3il-Xo)S 1 ^^^^ 

ri2 f 13 r23 ' 

where Aq is a free numerical parameter. Note that at S = 1/3, eq. ([T5|) reduces to eq. ([T0| . as required. We have used 
Aq = 0.5 throughout our computations since this value gave the best results for the H~ system, but other values of Aq 
may result in better convergence of the perturbation series for Ps^ . The effective potential, eq. (|14p , is minimized at 

fi2 = 0.773603828324(1), 
fi3 = r23 = 0.5866862922582(4), (16) 

with the minimum value of Vb = —1.381325607963162(1). The errors are estimated by performing the calculation 
again with higher precision and smaller tolerances. Convergence is ensured by restarting the minimization from a 
slightly perturbed location. We note that this result corresponds to the energies rescaled by 6^. Thus, to compare 
with the physical value, we have to divide this result by 3^ = 9, obtaining the first estimate of the binding energy 
~ —0.15 a.u., to be compared with the known value (see Table HB of about —0.26 a.u. 

The static 6 = limit is the zeroth order in 1/D expansion but, without the kinetic energy, it does not allow us 
to generate further orders in the perturbation expansion. In order to construct such an expansion, we consider the 
next simplest case, the harmonic approximation to the potential. This will yield a complete set of states that can be 
used to generate an expansion. The natural expansion parameter for eq. is 6^/^. This follows from the dominant 
balance argument applied to the Schrodinger equation. One finds that, for 5^/^, the harmonic terms in the expansion 
of the potential are of the same order as the constant coefficient terms in the kinetic energy expansion. 

Details of the procedure used to construct the expansion are described in the Appendix. The summation of the 
resulting series in powers oi 6 = 1 /D is complicated by the fact that the expansion is divergent at high orders due to 
a singularity at S = [s^ , so we expect the convergence of the naive summation 



fe=0 



to be slow. In the above expression Eq = Vq and E/- = £2/0-2 for fc > 0, where ek are expansion coefficients of the 
rescaled energy that appears in eq. ([T^ . There are also poles at _D = 1 that slow down the asymptotic convergence 
of the expansion at low values of D. A better estimate for E can be obtained by subtracting these poles from the 
expansion. To this end, the residues of the poles must be determined. Following |23l . Issj we define 



E{S) = 6^ 

where 



a_2 



Els'" 



(i-sy i-s ^ 

^ ' k=0 



(18) 



El = Ek-{k + l)a^2-a-i. (19) 

The residue of the second order pole, a_2, corresponds to the ground state energy in the D = 1 limit (more precisely 
a_2 = AEd^i). We have computed it employing again the method described in Section |lll this time with the rescaled 
charges of electrons and the positron so as to satisfy eq. (fTS]). We find 



a_2 = -1.102499999999999(1), (20) 

which again (see TableU]) resembles a rational number, indicating that there are likely analytical solutions of the D = 1 
model even for an arbitrary charge of the positron (not necessarily equal in magnitude to that of the electron) . As in 
Table HI the uncertainty in this converged residue is due to finite precision, as was checked by using larger quadrature 
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To find the residue of the single pole, a_i, we subtract the double pole from both sides of eq. (|T8)) and multiply by 
1 — S. We get the condition 

oo 

a_i = lim y^{Ek - Ek^i - a^2)S'', (21) 

k=0 

where E^i — 0. In practice we only have a finite number of terms in the sums in equations ((T8|) and (|21l) . Pade 
approximants have been shown to work well for summing up 1 /D expansions (22l . [23l | . Using this method to compute 
the limit in eq. (f2T1) we get for Ps^ , 

a_i = 0.427(2). (22) 

This result was obtained with the first 21 terms in the sum in eq. (|2ip . The uncertainty in the coinputed value was 
estimated by varying the order of the Pade approximant for a_i as [N/M] [{N — 1)/{M + 1)] [3g|. If the result 
has converged, the order of the approximant should not matter (barring the introduction of spurious poles in the 
denominator of the approximant). We use this method to estimate the error for all quantities computed using Pade 
approximants. In our calculations we use the full unrounded result for a_i which gives a slightly worse result for the 
bound state energy than eq. (j22p . As has been noted in Ref. ^2^, this way of determining a_i is not very accurate. 
An exact value for a_i (in principle obtainable from expansions about D — 1) would improve the convergence of 
the 1/D expansion. For He we used a_2 = —3.15546 [3l[ to get a__i = 0.313(1) using an identical calculation with 
21 energy expansion coefficients. We then evaluated eq. ([TSl) (with the summation truncated again at 21 terms) at 
S — 1/3. For helium, this yields a ground state energy that agrees with the variational calculation of [l^l to five 
digits, which is consistent with the result of [2l| for this summation method and perturbation expansion cutoff. The 
same calculation for the positronium ion yields a five digit agreement with the results in [lol These results 

are summarized in Table [III 





Known energy (la. u. — ma^) 


1/D Expansion 


He 


-2.9037243770341196 [37] 


-2.90374(1) 




-0.2620050702329801 [llj 


-0.262005(2) 



Table II. Results of summation of the 1/D expansion using Fade summation with first and second order poles removed. The 
first 21 non-zero terms were used in the summation. 



Figure[3]shows the improvements to the energy that are obtained by summing more terms. Higher orders yield better 
accuracy despite the poor behaviour of the \/D expansion coefhcients (see Table Hill . In fact, the pole subtraction 
and Pade resummation described above are necessary to get a sensible answer. 

Aside from computing higher orders in perturbation theory, precision of the result may be improved by using a 
different summation method. For example, Ref. [2l| found that Pade-Borel summation gives better results for helium 
than Pade summation. 



IV. CONCLUSIONS 

We have investigated the viability of dimensional scaling for making accurate predictions for the positronium ion 
system. Equal masses and correlation strengths make Ps~ a good candidate for the dimensional scaling treatment. 
We considered the D = 1 limit and found that the Schrddinger equation can be reduced to a one dimensional integral 
equation. The numerical solution for the energy eigenvalue approaches a simple rational number suggesting the 
possibility of a completely analytical solution. While this energy is not physically relevant by itself, it can be used to 
accelerate the convergence of the 1/D perturbation series. 

We constructed such a perturbative series by expanding the solution of the full Schrddinger equation about the 
Z? — >■ oo limit. Each coefficient was computed exactly in the harmonic basis. To obtain an accuracy of five significant 
figures required expanding up to order 41 in perturbation theory. While the accuracy of the energy expansion at this 
order is not yet competitive with variational calculations, the present method provides a valuable alternative approach 
to few body systems. It can be used to check a variety of matrix elements that have previously been computed only 
variationally. 

In the future, higher orders in the 1/D perturbation series can be determined without sacrificing speed if the 
analytical expansions can be replaced with numerical evaluations of series coefficients through finite differencing. It 



7 




I 1 > > > 1 

5 10 15 20 25 

Order n 

Figure 3. (Color online) Number of accurate digits in the ground state energy, defined as — logj^Q[(i5 — -Bexact)/-E'exact] ^ ^ 
function of the number of terms in the summation of eq. (|18|) . 



would also be very valuable to establish how the convergence of this expansion depends on the value of the parameter 
Ao introduced in eq. ((T5|) . Finally, the accuracy of the obtained wave function should be determined by evaluating 
matrix element of various operators and comparing them with the variational approach. 
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Appendix A: Perturbative expansion in 1/D 

In this Appendix we describe how the coefficients of the 1/D expansion were determined. Our procedure follows 
the matrix method of ref . . In terms of the displacement coordinates Xi defined by 

ri2 = + S^/'^xi 
ri3 = + S^/'^X2 

r23--r23 + S^^^X3, (Al) 
(fij are the coordinates of the minimum of the effective potential, eq. (|16p ). the Schrodinger equation takes the form 

[ST + Feff )</) = (A2) 
The Hamiltonian is expanded in powers of (5^/^ such that 

oo 

^ = (A3) 

1=0 
oo 

1=0 
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Since the expansion is about the minimum of Kff , there is no hnear term in its expansion and Vi = 0. Also the kinetic 
energy T starts contributing only in the second order in (5^/^, so the energy e has the form 

oo 

e = Fo + ^^e.(5*/^ (A5) 

For second order in S^^"^ in eq. (jA2l) we have 



d_Y f_d_Y f_d_^'' 
2 ..^ /'''^ KdxJ [3x2) \dx3, 



To--- ^ Ujk. 



i+j+k=2 
i+j+k=2 

where Ujk and Vijk are expansion coefficients that are functions of fij . This order in perturbation theory corresponds to 
three coupled harmonic oscillators. To solve the Schrodinger equation they need to be decoupled. This procedure yields 
normal mode frequencies uji and the corresponding normal coordinates g^, related to Xi by a linear transformation S, 

3 

q. = Y,iS-%^j- (A7) 



In terms of coordinates qi , 



1 A 92 



1 ^ 

1^2 = i'ooo + ^^w-g-. (A9) 



2 

i=l 



Defining 



the Hamiltonian can be written as 



H,^T, + V,+2, (AlO) 



H = ST + Vcff^Vo + HiS''^ . (All) 



i=0 

Next we consider the wave function expansion 



Y.4>^5''^■ (A12) 



i=0 

Without loss of generality (f>i can be normalized as 

{U^j) = 5^^,. (A13) 

Collecting like powers of 5^ in eq. (|A2p yields 

p 



^(il,-e,)<^p-, =0. (A14) 



1=0 

The p = Q order equation is a system of three independent harmonic oscillators with the solution 

'/'o = V (91)^1^2 (92)^1/3 (93), (A15) 

£0 = wooo + ^ Y^^^ ^" (A16) 
i=i ^ ^ 
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with 



K{qi) 



V ^ \r2yv\ 



(A17) 



where is the i^'th Hermite polynomial. For the ground state, Vi = 0. 

To compute further orders in the perturbation expansion, </)j from eq. (jA12|) are projected onto the harmonic 
oscillator basis 



(A18) 



Here jo^^'^^'^^ are the expansion coefficients. The advantage of using the Hermite function basis is that only a finite basis 
at every order of perturbation theory is needed, since the perturbations are polynomials in g^. Thus the perturbation 
expansion coefficients can be computed exactly. We note that oa*!*^'*^ = 5o,n'5o,i2<5o,i3- Equation (|Af3|) then implies 
that for any p > 



,000 



0. 



(A19) 



The matrix elements of the operators Hj defined in eq. (|Af 0[) are computed by noting that each Hj is a sum of terms 
of the form 



lo l'^ 

<?1 92 93 



_d_ 

dqi 



_d_ 

dq2 



d_ 



(A20) 



where ai + a2 + Q^s = 2 for the kinetic terms and ai = for terms co ming from VcQ. The matrix elements of qi and 
_d_ - - - - ' ^ 



are derived from the recurrence relations of the Hermite functions [20|, 



% = 



_d_ 
dqi 



/ VT \ 

\/T V2 
V2 V3 ••• 




V ■■) 

/ \/T 

-y/l \/2 

~V2 V3 ••■ 
-V3 



(A21) 



(A22) 



so Hj is a linear combination of direct products of such matrices. We denote the matrix representation of Hj by 
Hj , and by a.j the tensor with elements j-a*!*^*^ in the harmonic basis. Finally we derive the recursion relations for 
computation of the energy and wave function expansion coefhcients. First, we rewrite eq. (jAf4[) in the harmonic basis 



^(Hj - ei)ap_i = 0, 



(A23) 



i=0 



and then contract with ao and solve for e„, which yields 



= ao ^ HiBp- 



(A24) 



To compute the wave function expansion coefficients we need the pseudo- inverse K of the operator Hq — eo, defined 
component-wise as 



K 



fclfe2fe3 







iff^i^^^^'s^:^ otherwise 



(A25) 
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p 


Ep for Ps 


Cp for He 





-1.185438078904337(1) EO 


-2.423036748379509(1) El 


2 


-2.78770519314798(1) EO 


-3.544873487874171(2) El 


4 


-7.2695509791874(1) EO 


-5.56025516084019(2) El 


6 


6.7347904088005(2) El 


-2.174685942637(1) El 


8 


-2.1412953632562(1) E3 


-3.30958097736(1) E2 


10 


7.8884951280128(7) E4 


5.2508188805(2) E2 


12 


-3.519438146299(1) E6 


4.0504015254(2) E4 


14 


1.842029153744(1) E8 


-1.7333557830(1) E6 


16 


-1.107117203726(1) ElO 


5.6857880174(1) E7 


18 


7.51651405108(1) Ell 


-1.77525788344(2) E9 


20 


-5.69017274005(1) E13 


5.528541546(1) ElO 


22 


4.75290730575(1) E15 


-1.732045588(2) E12 


24 


-4.34262756760(1) E17 


5.409411228(4) E13 


26 


4.30867078063(3) E19 


-1.638399579(2) E15 


28 


-4.6135902290(1) E21 


4.4914966(2) E16 


30 


5.3029755999(1) E23 


-8.6523653(2) E17 


32 


-6.512769125044(1) E25 


-1.289717(2) E19 


34 


8.5114570184(2) E27 


3.152243(3) E21 


36 


-1.17941098599(2) E30 


-2.796124(6) E23 


38 


1.7272349225(1) E32 


2.031490(1) E25 


40 


-3.253630209(1) E34 


-1.705787919969(3) E33 



Table III. 1/D energy expansion coefficients in eq. (|A5|I . in units of ma^ . Terms with odd p vanish. Letter E indicates powers of 
10 multiplying the entries. The uncertainty in each coefficient was estimated using a similar calculation with 20 digit precision 
as described in the text. 



The operator K is defined such that K(Ho — eo) = 1 everywhere except for the subspace spanned by the harmonic 
ground state wave function (/)o, where the inverse of Hp — eo would be undefined and it is convenient to choose K =0. 
Contracting K with eq. (jA23|) gives 

p 

ap = K^(e, -H,)ap_,. (A26) 

i=l 

Together equations ()A24p and (IA26P allow us to compute the ground state energy and wave function to any order. 

We implemented the steps required to compute the 1/D expansion to arbitrary order in Mathematica [ssj and in 
CH — h. The determination of the Taylor expansion coefficients of the Hamiltonian is done with Mathematica. The 
computation of the perturbation series (eqs. (jA24[) and (|A26p ) is done in C++, for its speed of operations with large 
arrays (corresponding to the various tensor contractions in these equations). We have computed 20 1/D expansion 
coefficients (which required expanding up to order 41 in perturbation theory). The result is presented in Table Hill 
Also in this table are the corresponding coefficients for helium (from an identical calculation) , which agree to at least 
five significant figures with Table I of [2l| (after accounting for a difference in units, which amounts to diving by 
Z'^ = 4) and serve as a check of our calculations. Note that the coefficients in Table HlTl become large at high orders. 
This is due to the essential singularity at 6 = 0. The nature of this singularity has been investigated in Ref. (2]| . 

Due to the large number of algebraic operations required to generate the 1 /D expansion we need to check for round 
off error in our coefficients; one way to do this is to repeat the C++ computation at higher precision (there should 
be no need to redo the Mathematica part, since Mathematica does arbitrary precision computations by default, as 
long as one does not invoke numerical solvers). We have implemented a version of the CH — h code using the arbitrary 
precision arithmetic package ARPREC [s^ . The relative effect of round off error is shown in Fig. |4l We see that 
the error introduced by finite precision arithmetic is much smaller than the accuracy of the final ground state energy 
obtained by resumming the 1 /D expansion. Higher order calculations will require better precision when the fractional 
error becomes of the same order as the accuracy required. 
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Order n 

Figure 4. (Color online) Fractional error in the 1/D expansion coefficients defined as - E^^^^) / E^^°''\ as a function of 

the order n. The coefficients En^'^ were obtained using the standard double precision arithmetic (~ 16 digits of precision), 
while E^'^^ were obtained using the ARPREC arbitrary precision library (with 20 digits of precision). 
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